function ap_measure(varargin)
% AP_measure - measure action potential shape with "p/4" protocol subtraction
% 4/19/2007 P. Manis
%
% based on Analysis.m - plugin analysis template, Jan 2007
%
% Analysis is a template for data analysis routines to be used in the datac
% analysis program. This template shows how to write an analysis function,
% and how to allow the function to do an analysis on records selected in
% the main database window.
%%
% The template can be used as a starting point. If you change this, first
% change the name of the function as listed above to something reasonable.
% You do not have to change the name of the local function 'analysis2'.
% This code is more extensively documented than usual, to help you
% understand how it works. When you make your own routines, you may wish to
% remove some of these comments.
%
% Paul B. Manis, Ph.D. pmanis@med.unc.edu
%
% Variable and function descriptions:
% nargin: nargin is the count of the number of arguments that appeared
% in the function call to 'analysis'. nargin is provided by matlab when we
% call the routine. We use the 'varargin' (variable
% argument input list) for flexibility in writing the analysis routine so
% that we can call it from the command line for testing, or use it as a
% "plugin" with no arguments to do automated analysis.
%
% varargin: This is the matlab name for a variable argument list. A
% variable argument list is a cell array. Cell arrays can hold any data
% type, so you can mix text, numbers and structures (or other cell arrays)
% in the call. The array is indexed according to the position of the
% argument on the calling command line.
%
% test_mode: this is a flag that will be checked in the local routine. If
% it is 0, we access the data in the database for the analysis. If it is
% not, we use it as a directive to generate test data to verify the
% function of the routine.
%
% getmainselection: This is a function that we provide outside this
% function, to return the indices of the selected elements of the database
% (e.g., the lines in the main screen of datac when the database is shown).
% The return is a one-dimensional numeric array. Typically, I call this
% array 'sf' for 'selection'.
%
% getplotflag: This is a function that we provide outside this function. It
% reads the status of the checkbox on the datac window labeled "plot".
% You can use this flag to control whether the data is plotted on the
% display, or more frequently, whether it is sent automatically to the
% default printer during analysis. This allows you to run the analysis for
% test purposes without printing it, and then to rerun it with the printer
% enabled when everything is working so you can have hardcopy.
%
% analysis2: this is the routine that does the work. Functions called in
% 'analysis2' are described therein.

if(nargin == 0) % check for arguments in the call
    test_mode = 0; % if none, provide default values
else
    test_mode = varargin{1}; % arguments are in a cell array - extract them
end;
sf = getmainselection; % call the routine that lists the selected record indices
if(sf > 0) % make sure there's a selection
    pflag = getplotflag; % also get the plot flag checkbox status
    QueMessage('AP_measure', 1); % clear the que and show a message
    for i = 1:length(sf) % loop through the selected database entries
        analysis2(sf(i), pflag, test_mode); % call the analysis routine with argumetns
    end;
end;

return;


%------------------------------------------------------

function analysis2(sf, plot_flag, testmode)
% analysis2
%
% This routine does the actual work.
% The global variables below are used to acces the data and data structures
% in datac. ALLCH is a cell array matrix that holds the data for all channels.
% You really only need to access it if you have more than two channels, as
% same data for the first two channels is held in the VOLTAGE and CURRENT
% arrays (under the control of the mode flag - whether current or voltage
% clamp).
% DFILE is a structure that contains information about the most recently
% read data records. Some of the information is reduntantly stored in the
% CONTROL structure also.
% CONTROL is the central database of the whole datac program. CONTROL is an
% sturctured array. Each element of the array corresponds to an entry in
% the database - usually a group of records or a "block" of data from a
% cell. The elements of CONTROL contain information about the file for the
% data, the records, the gains, analysis parameters, and it can contain
% result information. The CONTROL structure is saved as the database. You
% can add fields to the CONTROL structure to store unique results.
%
% functions called:
%
% analysis_setup: This routine accesses the file pointed to by the selected
% (sf) record. It compares the current DFILE information with the database
% information in CONTROL(sf), and decides whether it needs to load a
% different set of data into the ALLCH, VOLTAGE and CURRENT arrays. If it
% does, it attempts to do so. If it is successful (finds the data and reads
% it without error), the return variable err will be 0, and DFILE will
% contain the new information about the data. If it fails for any reason
% (the data is not found at the path identified, or the file is corrupt),
% then err will be 1 and the program will return. DFILE will not be
% changed.
%
% analysis_setup2: This routine does common parsing of paramteres from the
% data and the stimulus structure. This parsing used to be done in-line,
% but it is much cleaner to just include it here. Note that the returned
% variables are just listed out.
%
% QueMessage: this routine prints a message to the mesage window of datac.

%
%

%
global   VOLTAGE 
global DFILE  CONTROL

[DFILE, err] = analysis_setup(DFILE, sf);
if(err ~= 0)
    return;
end;

[protocol, rate, records, pts, frec, lrec, time] = analysis_setup2(DFILE, sf);
QueMessage('AP_meaure Setup done');

% end of required code to read in the data and get parameters.
if(testmode)
end;


% min_isi = 10;
% 
% [stim, df] = block_info(sf);

do_spike = 1;

a=VOLTAGE(1:floor(records/2),:);

b=mean(VOLTAGE(floor(records/2):records,:));
rmp = mean(b(1:10));
CONTROL(sf).Rmp = rmp; % set RMP here - it is not set for this block elsewhere
b=b-mean(rmp); % zero baseline so we don't upset the rmp measure
b=b*4;
a4 = zeros(size(a));
for i = 1:size(a,1)
    a4(i,:) = a(i,:) +b; % subtraction of passive from active
end;
tsw1 = 10;
tsw2 = 0.75*max(max(time));
spike_thresh = 0;

QueMessage('AP_measure - Finding spikes');

[first_spike, first_isi, nr_spikes, spike_train]=find_spikes2(DFILE, a4, tsw1, tsw2, spike_thresh);

v = size(spike_train);
ISILAT(1:v(2))=deal(struct('nspikes', 0, 'isi', 0, 'lat', 0, 'source' , 0));
% initialize the ISILAT structure
for i=1:v(2)
    ISILAT(i).nspikes = length(spike_train(i).latency);
    ISILAT(i).isi= diff(spike_train(i).latency);
    ISILAT(i).lat = spike_train(i).latency;
    ISILAT(i).source = spike_train(i).source; % save the source trace number too
end

% deadwin = CONTROL(sf).deadwin; % dead window for most analyses EXCEPT spikes, which have their own values.

% analyze the spikes - hw, dv/dt, amplitude from rest, etc.
QueMessage('AP_mesure - Analyzing Spikes');
% initialize the structure
spike.timestamp = datestr(now);
spike.edit = 0; % not editable...
spike.n = 0;
spike.measthresh = NaN; % mean(spikesumm.measthresh); % fraction of height from rest to spike peak where measurement was made
spike.vthresh = NaN; % mean(spikesumm.vthr); % true threshold (2nd derivative method)
spike.vthrmin = NaN; % min(spikesumm.vthr); % min value (first spike...)
spike.deadwin = NaN; % mean(spikesumm.deadwin);
spike.minisi = NaN; % mean(spikesumm.minisi);
spike.max= NaN; % mean(spikesumm.vmax); spike.max_sd = std(spikesumm.vmax);
spike.apht = NaN; % mean(spikesumm.apht); spike.apht_sd = std(spikesumm.apht);
spike.half_ht = NaN; % mean(spikesumm.half_ht); spike.half_ht_sd = std(spikesumm.half_ht);
spike.min=NaN; % mean(spikesumm.vmin); spike.min_sd = std(spikesumm.vmin);
spike.dvmax = NaN; % mean(spikesumm.dvmax); spike.dvmax_sd = std(spikesumm.dvmax);
spike.dvmin = NaN; % mean(spikesumm.dvmin); spike.dvmin_sd = std(spikesumm.dvmin);
%      dvvmax = spikesumm.dvvmax(find(~isnan(spikesumm.dvvmax)));
spike.dvvmax = NaN; % mean(dvvmax); spike.dvvmax_sd = std(dvvmax);
%      dvvmin = spikesumm.dvvmin(find(~isnan(spikesumm.dvvmin)));
spike.dvvmin = NaN; % mean(dvvmin); spike.dvvmin_sd = std(dvvmin);
%     hwa = spikesumm.hwa(find(~isnan(spikesumm.hwa))); hwb = spikesumm.hwb(find(~isnan(spikesumm.hwb)));
spike.hwa = NaN; % mean(hwa); spike.hwa_sd = std(hwa); % eliminate the ones where coulnd't get a measurement
spike.hwb = NaN; % mean(hwb); spike.hwb_sd = std(hwb);
spike.ahpdepth =  NaN; % mean(spikesumm.ahpdepth); spike.ahpdepth_sd = std(spikesumm.ahpdepth);

CONTROL(sf).spike=spike; % set up the initial array

nspike = size(ISILAT);
min_isi = 10;
if(nspike(2) > 0 && (do_spike == 1))
    spikesumm = spike_sum(DFILE, a4, ISILAT, tsw1, CONTROL(sf).Rmp, ...
        CONTROL(sf).filename, CONTROL(sf).recbeg, CONTROL(sf).spike.deadwin, ...
        min_isi, tsw1, tsw2, plot_flag);
    % summary spike information:
    if(~isempty(spikesumm))
        n = size(spikesumm.vmax);
        spike.timestamp = datestr(now);
        spike.edit = 0; % not editable...
        spike.n = n(2);
        spike.measthresh = mean(spikesumm.measthresh); % fraction of height from rest to spike peak where measurement was made
        athr = find(~isnan(spikesumm.vthr));
        if(isempty(athr))
            spike.vthresh=NaN;
            spike.vthrmin=NaN;
        else
            spike.vthresh = mean(spikesumm.vthr(athr)); % true threshold (2nd derivative method)
            spike.vthrmin = min(spikesumm.vthr(athr)); % min value (usually the first spike...)
        end;
        spike.deadwin = mean(spikesumm.deadwin);
        spike.minisi = mean(spikesumm.minisi);
        spike.max=mean(spikesumm.vmax); spike.max_sd = std(spikesumm.vmax);
        spike.apht = mean(spikesumm.apht); spike.apht_sd = std(spikesumm.apht);
        spike.half_ht = mean(spikesumm.half_ht); spike.half_ht_sd = std(spikesumm.half_ht);
        spike.min=mean(spikesumm.vmin); spike.min_sd = std(spikesumm.vmin);
        spike.dvmax = mean(spikesumm.dvmax); spike.dvmax_sd = std(spikesumm.dvmax);
        spike.dvmin = mean(spikesumm.dvmin); spike.dvmin_sd = std(spikesumm.dvmin);
        dvvmax = spikesumm.dvvmax(find(~isnan(spikesumm.dvvmax))); %#ok<FNDSB>
        spike.dvvmax = mean(dvvmax); spike.dvvmax_sd = std(dvvmax);
        dvvmin = spikesumm.dvvmin(find(~isnan(spikesumm.dvvmin))); %#ok<FNDSB>
        spike.dvvmin = mean(dvvmin); spike.dvvmin_sd = std(dvvmin);
        hwa = spikesumm.hwa(find(~isnan(spikesumm.hwa))); hwb = spikesumm.hwb(find(~isnan(spikesumm.hwb))); %#ok<FNDSB>
        spike.hwa = mean(hwa); spike.hwa_sd = std(hwa); % eliminate the ones where coulnd't get a measurement
        spike.hwb = mean(hwb); spike.hwb_sd = std(hwb);
        spike.ahpdepth =  mean(spikesumm.ahpdepth); spike.ahpdepth_sd = std(spikesumm.ahpdepth);
        CONTROL(sf).spike=spike; % make it a substructure
        % spike train stuff...
        if(plot_flag >= 0)
            h = findobj('Tag', 'CC_IV_spiketrain'); % check for pre-existing window
            if(isempty(h)) % if none, make one
                h = figure('Tag', 'CC_IV_spiketrain', 'Name', 'spike train Analysis', 'NumberTitle', 'off');
            end
            figure(h); % otherwise, select it
            clf; % always clear the window...
            subplot(2,2,1);
            plot(spikesumm.lat, spikesumm.dvmax, 'gx');
            xlabel('T (ms)');
            ylabel('dvmax (mV/ms) - rise');
            subplot(2,2,2);
            plot(spikesumm.lat, spikesumm.dvmin, 'rx');
            xlabel('T (ms)');
            ylabel('dvmin (mV/ms) - fall');

            subplot(2,2,3);
            plot(spikesumm.lat(find(~isnan(spikesumm.hwa))), hwa, 'go'); %#ok<FNDSB>
            xlabel('T (ms)');
            ylabel('hwa (ms) - rising half-width');
            subplot(2,2,4);
            plot(spikesumm.lat(find(~isnan(spikesumm.hwb))), hwb, 'ro'); %#ok<FNDSB>
            xlabel('T (ms)');
            ylabel('hwb (ms) - falling half-width');

        end;
    end;
end;
CONTROL(sf).spike=spike;

